Method and Apparatus for Efficient Three-Dimensional Contouring of Medical Images

ABSTRACT

A technique is disclosed for generating a new contour and/or a 3D surface such as a variational implicit surface from contour data. In one embodiment, a point reduction operation is performed on data sets corresponding to any combination of transverse, sagittal, or coronal contour data prior to processing those data sets to generate a 3D surface such as a variational implicit surface. A new contour can also be generated by the intersection of this surface with an appropriately placed and oriented plane. In this manner, the computation of the variational implicit surface becomes sufficiently efficient to make its use for new contour generation practical.

CROSS-REFERENCE AND PRIORITY CLAIM TO RELATED PATENT APPLICATIONS

This application is a divisional of U.S. patent application Ser. No.11/848,624 filed Aug. 31, 2007, published as U.S. patent applicationPub. 2009/0060299, now U.S. Pat. No. ______, the entire disclosure ofwhich is incorporated herein by reference.

This application is also related to U.S. patent application Ser. No.______ filed this same day, and identified by Thompson Coburn AttorneyDocket Number 37340-100047, which is a divisional U.S. patentapplication Ser. No. 11/848,624 filed Aug. 31, 2007, published as U.S.patent application Pub. 2009/0060299, now U.S. Pat. No. ______.

FIELD OF THE INVENTION

The present invention pertains generally to the field of processingmedical images, particularly generating contours for three-dimensional(3D) medical imagery.

BACKGROUND AND SUMMARY OF THE INVENTION

Contouring is an important part of radiation therapy planning (RTP),wherein treatment plans are custom-designed for each patient's anatomy.Contours are often obtained in response to user input, wherein a usertraces the object boundary on the image using a computer workstation'smouse and screen cursor. However, it should also be noted that contourscan also be obtained via automated processes such as auto-thresholdingprograms and/or auto-segmentation programs.

FIG. 1 depicts an exemplary GUI 100 through which a user can view andmanipulate medical images. The GUI 100 includes frame 102 correspondingto the transverse (T) viewing plane, frame 104 corresponding to thecoronal (C) viewing plane, and frame 106 corresponding to the sagittal(S) viewing plane. Within frame 102, an image slice of a patient thatresides in a T plane can be viewed. Within frame 104, an image slice ofa patient that resides in a C plane can be viewed. Within frame 106, animage slice of a patient that resides in an S plane can be viewed. Usingwell-known techniques, users can navigate from slice-to-slice andviewing plane-to-viewing plane within GUI 100 for a given set of imageslices. It can also be noted that the upper right hand frame of GUI 100depicts a 3D graphics rendering of the contoured objects.

FIG. 2( a) illustrates an exemplary patient coordinate system withrespect to a radiotherapy treatment machine that is consistent with thepatient coordinate system defined by the IEC 61217 Standard forRadiotherapy Equipment. As can be seen, the patient coordinate system isa right-hand coordinate system such that if a supine patient is lying ona treatment couch with his/her head toward the gantry, the positivex-axis points in the direction of the patient's left side, the positivey-axis points in the direction of the patient's head, and the positivez-axis points straight up from the patient's belly. The origin of thiscoordinate system can be offset to the origin of the image data understudy.

FIG. 2( b) defines the T/S/C viewing planes with respect to the patientcoordinate system of FIG. 2( a). As is understood, a plane in the Tviewing plane (the xz-viewing plane) will have a constant value for y, aplane in the S viewing plane (the yz-viewing plane) will have a constantvalue for x, and a plane in the C viewing plane (the xy-viewing plane)will have a constant value for z.

Returning to the example of FIG. 1, the image data within GUI 100depicts a patient's prostate 110, bladder 112, and rectum 114. Asindicated above, an important part of RTP is the accurate contouring ofregions of interest such as these.

Current RTP software typically limits contour drawing by the userthrough GUI 100 to T views (views which are perpendicular to thepatient's long axis) as the T images usually have the highest spatialresolution, the T images are the standard representation of anatomy inthe medical literature, and the T contours are presently the only formatdefined in the DICOM standard. The two other canonical views—the S and Cviews—can then be reconstructed from the columns and rows, respectively,of the T images.

When generating 3D surfaces from image slices, conventional softwareprograms known to the inventor herein allow the user to define multipleT contours for a region of interest within an image for a plurality ofdifferent T image slices. Thereafter, the software program is used tolinearly interpolate through the different T contours to generate a 3Dsurface for the region of interest. However, the inventor herein notesthat it is often the case that a plane other than a T plane (e.g.,planes within the S and/or C viewing planes) will often more clearlydepict the region of interest than does the T plane. Therefore, theinventor herein believes there is a need in the art for a robust 3Dcontouring algorithm that allows the user to define input contours inany viewing plane (including S and C viewing planes) to generate a 3Dsurface for a region of interest and/or generate a new contour for theregion of interest.

Further still, the inventor herein believes that conventional 3D surfacegeneration techniques, particularly techniques for generatingvariational implicit surfaces, require unacceptably long computationaltimes. As such, the inventor herein believes that a need exists in theart for a more efficient method to operate on contours in threedimensions.

Toward these ends, according to one aspect of an embodiment of theinvention, disclosed herein is a contouring technique that increases theefficiency of 3D contouring operations by reducing the number of datapoints needed to represent a contour prior to feeding those data pointsto a 3D contouring algorithm, wherein the 3D contouring algorithmoperates to generate a 3D surface such as a variational implicit surfaceor process the reduced data points to generate a new contour in a newplane via an interpolation technique such as B-spline interpolation. Thedata points that are retained for further processing are preferably aplurality of shape-salient points for the contour. In accordance withone embodiment, computed curvature values for the data points are usedas the criteria by which to judge which points are shape-salient. Inaccordance with another embodiment, computed scalar second derivativevalues are used as the criteria by which to judge which points areshape-salient. In accordance with yet another embodiment, the DeBoorequal energy theorem is used as the criteria by which to judge whichpoints are shape-salient.

According to another aspect of an embodiment of the invention, disclosedherein is a contouring technique that operates on a plurality of datapoints, wherein the data points define a plurality of contourscorresponding to a region of interest within a patient, each contourbeing defined by a plurality of the data points and having acorresponding plane, wherein the plurality of data points are reduced asdescribed above and processed to find the reduced data points thatintersect a new plane, and wherein B-spline interpolation is used tointerpolate through the points of intersection to generate a new contourin the new plane. This embodiment can operate on a plurality of contoursdrawn by a user in the S and/or C viewing planes to generate a T contourin a desired T plane. The point reduction operation performed prior tothe B-spline interpolation improves the efficiency of the B-splineinterpolation operation.

While various advantages and features of several embodiments of theinvention have been discussed above, a greater understanding of theinvention including a fuller description of its other advantages andfeatures may be attained by referring to the drawings and the detaileddescription of the preferred embodiment which follow.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 depicts an exemplary graphical user interface (GUI) for acontouring program, wherein the GUI displays 3D patient ComputedTomography (CT) data in separate planar views;

FIG. 2( a) depicts an exemplary patient coordinate system with respectto a radiotherapy treatment machine;

FIG. 2( b) depicts the T, S, and C view planes for the patientcoordinate system of FIG. 2( a);

FIG. 3 depicts an exemplary contour specified by cubic B-splines;

FIG. 4 depicts an exemplary contour and its approximation using cubicB-spline interpolation using eight points on the original contour;

FIGS. 5( a) and (b) depict exemplary computing environments on whichembodiments of the present invention can be realized;

FIG. 6 depicts an exemplary process flow for generating a new contourfrom a plurality of input points that represent contours;

FIG. 7 depicts a plot of curvature versus arc length for sampled pointsfrom the exemplary contour of FIG. 4 for two finite difference intervalsettings;

FIG. 8 depicts the exemplary contour of FIG. 4 reconstructed usingB-spline interpolation for several different finite difference intervalsettings with respect to the curvature computation;

FIG. 9 depicts plots of scalar second derivative versus arc length fortwo settings of the finite difference interval for the exemplary contourof FIG. 4;

FIG. 10 depicts the exemplary contour of FIG. 4 reconstructed usingB-spline interpolation for several different finite difference intervalsettings with respect to the scalar second derivative computation;

FIG. 11 depict a plot of the DeBoor energy versus arc length and a plotof the cumulative distribution of energy for the exemplary contour ofFIG. 4;

FIG. 12 depicts three DeBoor energy reconstructions of the exemplarycontour of FIG. 4;

FIGS. 13( a)-(c) depict a graphical reconstruction of T contours usingB-spline interpolation in accordance with an embodiment of theinvention;

FIG. 14 depicts an exemplary process flow for generating a variationalimplicit surface from a plurality of input points that representcontours;

FIG. 15 depicts, for a pair of actual manually drawn contours, thenon-uniform data sampling, the results of the second derivative shapeanalysis, and the construction of a complete set of constraints forinput to an implicit function computation;

FIG. 16 depicts a variational implicit surface produced from thecontours of FIG. 15;

FIGS. 17-19( b) demonstrate the application of an exemplary variationalimplicit surface method in the contouring of the prostate, bladder, andrectum shown in FIG. 1;

FIG. 20 depicts an exemplary process flow for operating on input pointsets wherein one or more of the input sets contains a small number ofcontour data points; and

FIG. 21 depicts an exemplary process flow wherein both B-splinereconstruction of T contours and variational implicit surface generationfrom reduced data sets is employed.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

I. Contours:

The embodiments of the present invention address contours. Contours areplanar, closed curves C(x,y,z) which can be realized as sets ofnon-uniformly sampled points along the user-input stroke, {c₁, . . .,c_(M)} (or sets of points generated by an auto-thresholding and/orauto-segmentation program), wherein the individual points arerepresented by c_(i)=C(x_(i),y_(i),z_(i)), and wherein M is the numberof points in the contour. Points c_(i) in the T planes (xz-planes) havey constant, S contours (yz-planes) have x constant, and C contours(xy-planes) have z constant.

Contours can also be parameterized by a curve length u where the curve Cof length L is represented as C(x,y,z)=C(x(u),y(u),z(u))=C(u) where0≦u≦L and C(0)=C(L)

II. B-Spline Representation of Contours:

When contours exist as discrete points as noted above, it can be usefulto represent these points as samples on a continuous curve along whichone can interpolate the contour shape at any arbitrary point. B-splines,which can specify arbitrary curves with great exactness, can providesuch a representation for contours. (See Piegl, L. A., and Tiller, W.,The Nurbs Book, Springer, New York, 1996, the entire disclosure of whichis incorporated herein by reference). The B-spline description of acurve depends on (1) a set of predefined basis functions, (2) a set ofgeometric control points, and (3) a sequence of real numbers (knots)that specify how the basis functions and control points are composed todescribe the curve shape. Given this information, the shape of C(u) canbe computed at any u. Alternatively, given points u′ sampled along C(u),one can deduce a set of B-spline control points and corresponding knotsthat reconstruct the curve to arbitrary accuracy. Thus, B-splines can beused to interpolate curves or surfaces through geometric points or toapproximate regression curves through a set of data points.

B-splines form piecewise polynomial curves along u, delimited by theknots u_(i),i=0, . . . , m into intervals in which subsets of the basisfunctions and the control points define C(u). The m+1 knots U={u₀, . . .,u_(m)} are a non-decreasing sequence of real numbers such thatu_(i)≦u_(i+1), for all i.

The p-th degree B-spline basis function, N_(i,p)(u), defined for thei-th knot interval, defines the form of the interpolation. The zero-thorder function, N_(i,0)(u), is a step function and higher orders arelinear combinations of the lower order functions. The construction ofbasis functions by recursion is described in the above-referenced workby Piegl and Tiller. A preferred embodiment of the present inventiondescribed herein employs cubic (p=3) B-splines.

Basis function N_(i,p)(u) is nonzero on the half-open interval[u_(i),u_(i+p+1)), and for any interval [u_(i),u_(i+1)) at most (p+1) ofthe basis functions, N_(i−p,p)(u), . . . ,N_(ix,p)(u), are nonzero. Ap-th degree, open B-spline curve C(u) with end points u=a,b is definedby

$\begin{matrix}{{{C(u)} = {\sum\limits_{i = 0}^{n}{{N_{i,p}(u)}P_{i}}}},{a \leq u \leq b}} & (1)\end{matrix}$

where the P_(i) are the (n+1) control points, the N_(i,p)(u) are thebasis functions, and the knot vector U is defined

$\begin{matrix}{U = \left\{ {\underset{p + 1}{\underset{}{a,\ldots \mspace{14mu},a}},u_{p + 1},\ldots \mspace{14mu},u_{m - p - 1},\underset{p + 1}{\underset{}{b,\ldots \mspace{14mu},b}}} \right\}} & (2)\end{matrix}$

where a≦u_(p+1)≦u_(p+2), . . . ≦u_(m−p−1)≦b. This defines an unclosedcurve with multiple knots at the end values a=u₀, . . .,u_(p);b=u_(m−p), . . . ,u_(m). For a spline of degree p with m+1 knots,n+1 control points will be required to specify the shape; for all splinegeometries p,n,m are related as

m=n+p+1.  (3)

Closed curves with coincident start and end points and with C²continuity (continuous curve with continuous first and secondderivatives) throughout are defined with uniform knot vectors of theform U={u₀,u₁, . . . ,u_(m)} with n+1(=m−p) control points defined suchthat the first p control points P₀,P₁, . . . ,P_(p-1) are replicated asthe last p control points P_(n−p−1), . . . ,P_(n) which for the cubic(p=3) case means that P₀=P_(n-2), P₁=P_(n-1),P₂=P_(n). This means thatthere are actually n+1−p unique control points, and that the knots thatare actually visualizable on a closed curve are the set u_(p), u_(p+1),. . . , u_(m−p−1).

FIG. 3 shows an exemplary continuous closed curve 300 specified by cubic(p=3) B-splines. Curve 300 is defined by five unique control pointsP₀,P₁, . . . ,P₄ as shown at the vertices of polygon 302. To generate aclosed curve with C² continuity everywhere, p=3 of the control pointsare replicated, making n=7, and since m=n+p+1, then m+1=12 uniformlyspaced knots are required, given by the vector

U=(0,1,2,3,4,5,6,7,8,9,10,11)

For fixed p,n,U, the curve shape 300 can be changed by moving one ormore of the control points P. The locations of the knots are shown asdots on curve 300, wherein the knots u₃-u₇ uniquely span the curve,wherein knots u₀-u₂ coincide with knots u₅-u₇, and wherein knots u₈-u₁₁coincide with knots u₃-u₆. Thus, as with the control points that must beduplicated for cyclic B-spline curves, so too must some of the knots beduplicated.

III. B-Spline Interpolation of Points in Contours:

A useful application of B-splines is to interpolate a smooth curvethrough a series of isolated points that represent samples of a curve.Global interpolation can be used to determine a set of control pointsgiven all the data in the input curve. (See Chapter 9 of theabove-referenced work by Piegl and Tiller). Suppose one starts with aset of points {Q_(k)},k=0, . . . ,n on the actual curve, and the goal isto interpolate through these points with a p-degree B-spline curve.Assigning a parameter value ū_(k) to each Q_(k) and selecting anappropriate knot vector U={u₀, . . . ,u_(n)}, one can then set up the(n+1)×(n+1) system of linear equations

$\begin{matrix}{Q_{k} = {{C\left( {\overset{\_}{u}}_{k} \right)} = {\sum\limits_{i = 0}^{n}{{N_{i,p}\left( {\overset{\_}{u}}_{k} \right)}P_{i}}}}} & (4)\end{matrix}$

where the n+1 control points P_(i) are the unknowns. The system can bere-written as

Q=AP  (5)

where the Q,P are column vectors of the Q_(k) and P_(i), respectively,and where A is the matrix of basis functions. This (n+1)×(n+1) linearsystem can be solved for the unknown control points P_(i)

P=A⁻¹Q  (6)

by factoring A by LU decomposition instead of inverting matrix A. (SeePress, et al., Numerical Recipes in C, 2^(nd) Edition, CambridgeUniversity Press, 1992; Golub, G. H. and Van Loan, C. F., MatrixComputations, The Johns Hopkins University Press, Baltimore, 1996, theentire disclosures of both of which are incorporated herein byreference). A higher quality reconstruction—end points joined with C²continuity—can be obtained by restricting curves to cubic (p=3) type andby specifying endpoint first derivatives. Defining the endpoint tangentvectors D₀ at Q₀ and D_(n) at Q_(n), one constructs a linear system likeequation (5) but with two more variables to encode the tangentinformation resulting in a (n+3)×(n+3) system. The tangents are added tothe system with the equations

$\begin{matrix}{{P_{0} = {{Q_{0} - P_{0} + P_{1}} = {{{\frac{u_{4}}{3}D_{0}} - P_{n + 1} + P_{n + 2}} = {\frac{1 - u_{n + 2}}{3}D_{n}}}}}{P_{n + 2} = Q_{n}}} & (7)\end{matrix}$

that can be used to construct a tridiagonal system

$\begin{matrix}{\begin{bmatrix}{Q_{1} - {a_{1}P_{1}}} \\Q_{2} \\\vdots \\Q_{n - 2} \\{Q_{n - 1} - {c_{n - 1}P_{n + 1}}}\end{bmatrix} = {\begin{bmatrix}b_{1} & c_{1} & 0 & \ldots & 0 & 0 & 0 \\a_{2} & b_{2} & c_{2} & \ldots & 0 & 0 & 0 \\\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\0 & 0 & 0 & \ldots & a_{n - 2} & b_{n - 2} & c_{n - 2} \\0 & 0 & 0 & \ldots & 0 & a_{n - 1} & b_{n - 1}\end{bmatrix}\begin{bmatrix}P_{2} \\P_{3} \\\vdots \\P_{n - 1} \\P_{n}\end{bmatrix}}} & (8)\end{matrix}$

that can be solved by Gaussian elimination. (See Chapter 9.2.3 of theabove-referenced work by Piegl and Tiller).

To demonstrate the interpolation of points representing a putativecurve, the inventor has sampled points from closed curves with random,but known, shapes, and reconstructed the random curves measuring theaccuracy as the mean squared error of the reconstructed curve versus theoriginal. FIG. 4 shows an example of a contour 400 from which eightQ_(k) points were selected, and from which a set of control points 412was computed using equation (8). Contour 400 represents arandomly-generated shape, wherein this shape is approximated by contour406 based on cubic B-spline interpolation through eight points Q_(k) onthe original contour 400. Shown at right in FIG. 4 is a superposition ofthe approximated contour 406 over the original contour 400. Also shownat right in FIG. 4 is the geometric polygon representation 410 of thespline reconstruction, wherein the newly-determined control points arethe vertices 412 of polygon 410, and wherein the knots 408 on the newcurve 406 are connected by the line segments to the control points 412.In this example, there are 11 (which equals (n+1)) actual control points412, which include eight unique control points and the p replicates.Because m=n+p+1=10+3+1=14 , the knot vector requires m+1=15 elements forwhich one can define a uniform (equal knot intervals) knot vector,U=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14) , which is the usual knotconfiguration for closed curves.

IV. Embodiments of the Invention:

FIGS. 5( a) and 5(b) depict exemplary computing environments in whichembodiments of the present invention can be realized. Preferably, aprocessor 502 is configured to execute a software program to carry outthe three-dimensional contouring operations described herein. Such asoftware program can be stored as a set of instructions on anycomputer-readable medium for execution by the processor 502. Theprocessor 502 receives as inputs a plurality of data points 500, whereinthese input data points 500 are representative of a plurality ofcontours. These input points can be defined manually by a computer user(e.g., by dragging a mouse cursor over a desired shape to define pointsfor an input contour) or automatically by an auto-thresholding and/orauto-segmentation process, as would be understood by those havingordinary skill in the art.

In the embodiment of FIG. 5( a), the input points 500 are representativeof a plurality of contours that reside in at least one viewing plane.Preferably, these contours are a plurality of S contours, a plurality ofC contours, or some combination of at least one S contour and at leastone C contour. Furthermore, as described hereinafter, the softwareprogram executed by processor 502 is preferably configured to generateone or more output contours 504 from the input points 500, wherein theoutput contour(s) 504 reside in a viewing plane that is non-parallel tothe at least one viewing plane for the contours of the input points 500.Preferably, the output contour(s) 504 is/are T contour(s).

FIG. 6 depicts an exemplary process flow for the FIG. 5( a) embodiment.At step 602, the software program receives the plurality of input datapoints 500. These input points 500 can be grouped as a plurality ofdifferent initial sets of data points, wherein each initial point set isrepresentative of a different contour, the different contours existingin at least one viewing plane. As explained above, these contourspreferably comprise (1) a plurality of S contours, (2) a plurality of Ccontours, or (3) at least one S contour and at least one C contour.

One observation that can be made from FIG. 4 is that the reconstructedcontour 406 does not match the original contour 400 everywhere. The fitcould be improved by sampling more points on the original contour 400,and, in the limit of all available points on the original contour 400,the cubic B-spline reconstruction would be exact. However, the inventorherein notes that some points corresponding to the original contour 400ought to be more important to contour shape than others—e.g., points atu where the curvature is large should convey more shape information thanany other points. Thus, in the interest of increasing the efficiencywith which B-spline interpolation can be performed and with whichinterpolated contours can be represented through data, the inventorherein notes that B-spline interpolation need not be performed on all ofthe input points 500 for a given contour. Instead, the input points 500in a given set of input points can be processed to generate a reducedset of input points on which the B-spline interpolation will beperformed (step 604). Preferably, step 604 operates to reduce the numberof input points 500 for a given contour by generating a reduced set ofinput points, wherein the reduced set comprises a plurality ofshape-salient points. As used herein, points which are representative ofa contour are considered “shape-salient” by applying a filteringoperation based on a shape-indicative metric to those points, examplesof which are provided below.

The inventor herein discloses three techniques that can be used toreduce the data points 500 to a plurality of shape-salient points.

According to a first technique of point reduction for step 604, theshape-salient points for each initial input point set are determined asa function of computed curvature values for a contour defined by thepoints 500 within that initial input point set. The curvature isrepresentative of the speed at which curve C(u) changes direction withrespect to increasing u, wherein u represents the distance along curveC(u) beginning from an arbitrary starting point or origin. The curvatureof a plane curve is defined as:

$\begin{matrix}{{\kappa = \frac{{x^{\prime}y^{''}} - {y^{\prime}x^{''}}}{\left\lbrack {x^{\prime 2} + y^{\prime 2}} \right\rbrack^{\frac{3}{2}}}},} & (9)\end{matrix}$

where x′=dx/du, x″=d²x/du², etc. are derivatives computed by finitedifferences on uniform u−intervals along C(u), and where the x,y valuescorrespond to points which are representative of the input contour. (SeeDoCarmo, M., Differential Geometry of Curves and Surfaces, PrenticeHall, New York, 1976; Thomas, J. W., Numerical Partial DifferentialEquations—Finite Difference Methods, Springer, New York, 1995, theentire disclosures of which are incorporated herein by reference).

Preferably, step 604 takes points u* at peak values of κ(u)

u*=arg max_(u)κ(u)  (10)

These points u*, which contribute most importantly to the shape of acurve, are saved for reconstruction of the contour through B-splineinterpolation. It should be noted that because of the cyclic nature ofthe data in u (since 0≦u≦L and C(0)=C(L)), when computing the argmaxfunction over intervals u, one can let the intervals span the origin 0and then reset the computation for intervals placed at L+a to a or −a toL−a. To accomplish the use of uniform intervals u along C(u), one can(1) reconstruct each input contour via B-spline interpolation throughall of its raw input points, (2) step along the reconstructed contour inequal size steps that are smaller than the normal spacing among the rawinput points to generate the points which are fed to the curvaturecomputation of formula (9), and (3) apply the curvature computations offormulas (9) and (10) to thereby generate a set of reduced points fromthe original set of raw input points.

FIG. 7 shows plots of curvature versus arc length u for two settings ofthe finite difference interval for the random contour 400 of FIG. 4.Plot 700 shows the curvature versus arc length for finite differencesover 1.25% of the total curve length, and plot 702 shows the curvatureversus arc length for a finite differences of over 5.0% of the totalcurve length. By varying this interval size, one may detect more orfewer peaks along arc length u; the longer the interval, the moreapparent smoothing of the shape, thereby resulting in only the mostprominent peaks being detected. For example, in plot 702, only the mostprominent directional changes are apparent, leading to a selection oftwelve unique Q_(k) points for the B-spline interpolation analysis. Itshould also be noted that rather than using only maxima, step 604 canalso be configured to retain only those points for which the computedcurvature value exceeds a threshold value. As such, it can be seen thata variety of conditions can be used for determining how the curvaturevalues will be used to define the shape-salient points.

FIG. 8 shows the random shape 400 of FIG. 4 (shown with dotted lines inFIG. 8) reconstructed using B-spline interpolation for several intervalsettings. The number of points Q_(k) for each reconstruction is set bythe value of the finite distance interval and also depends on the shapecomplexity of the curve. For reconstructed contour 802, the number ofpoints Q_(k) was 12. For reconstructed contour 804, the number of pointsQ_(k) was 16. For reconstructed contour 806, the number of points Q_(k)was 29. As can be seen from FIG. 8, by using more points Q_(k), it ispossible to achieve arbitrarily high accuracy. However, it should alsobe noted that given the original number of samples for this example (atotal of 1,099 samples), the use of only 29 points represents asignificant compression in the number of points Q_(k) used for theB-spline interpolation while still providing good accuracy in thereconstruction. Optionally, user input can be used to define a settingfor the finite difference interval used in the curvature computations.Thus, the finite difference interval size in the curvature calculationcan serve as an adjustable tuning parameter which can be controlled todefine a quality and efficiency for the contour reconstruction.

According to a second technique of point reduction for step 604, theshape-salient points for each initial input point set are determined asa function of computed scalar second derivative a values (i.e., thescalar acceleration) for the motion of a point along C(u), which isdefined as

$\begin{matrix}{{a(u)} = \left\lbrack {\left( \frac{^{2}{x(u)}}{u^{2}} \right)^{2} + \left( \frac{^{2}{y(u)}}{u^{2}} \right)^{2}} \right\rbrack^{\frac{1}{2}}} & (11)\end{matrix}$

Preferably, step 604 takes points u* at peak values of a(u),

u*=arg max_(u)a(u)  (12)

Once again, the derivatives can be computed by finite differences onuniform u−intervals along C(u). Also, as noted above, because of thecyclic nature of the data in u (since 0≦u≦L and C(0)=C(L)), whencomputing the argmax function over intervals u, one can let theintervals span the origin 0 and then reset the computation for intervalsplaced at L+a to a or −a to L−a. As with the curvature calculationsdescribed above, to accomplish the use of uniform intervals u alongC(u), one can (1) reconstruct each input contour via B-splineinterpolation through all of its raw input points, (2) step along thereconstructed contour in equal size steps that are smaller than thenormal spacing among the raw input points to generate the points whichare fed to the scalar second derivative computation of formula (11) ,and (3) apply the scalar second derivative computation of formulas (11)and (12) to thereby generate a set of reduced points from the originalset of raw input points.

FIG. 9 shows plots of scalar second derivative versus arc length u fortwo settings of the finite difference interval for the random contour400 of FIG. 4. Plot 900 shows the scalar second derivative versus arclength for finite differences over 5.0% of the total curve length, andplot 902 shows the scalar second derivative versus arc length for afinite differences of over 20.0% of the total curve length. As with thecurvature method exemplified by FIG. 7, by varying the finite differenceinterval size, one may detect more or fewer peaks along arc length u;the longer the interval, the more apparent smoothing of the shape,thereby resulting in only the most prominent peaks being detected.Therefore, plot 902 would be expected to produce fewer Q_(k) points thanplot 900.

It should also be noted that rather than using only maxima, step 604 canalso be configured to retain only those points for which the computedscalar second derivative exceeds a threshold value. As such, it can beseen that a variety of conditions can be used for determining how thescalar second derivative values will be used to define the shape-salientpoints.

FIG. 10 shows the random shape 400 of FIG. 4 (shown with dotted lines inFIG. 10) reconstructed using B-spline interpolation for several intervalsettings. As with the reconstructions of FIG. 8, the number of pointsQ_(k) for each reconstruction is set by the value of the finite distanceinterval and also depends on the shape complexity of the curve. Forreconstructed contour 1002, the number of points Q_(k) was 10. Forreconstructed contour 1004, the number of points Q_(k) was 19. Forreconstructed contour 1006, the number of points Q_(k) was 44. As can beseen from FIG. 10, by using more points Q_(k), it is possible to achievearbitrarily high accuracy. Also, as with the reconstructions of FIG. 8,it should be noted that given the original number of samples for thisexample (a total of 1,099 samples), the use of only 44 points representsa significant compression in the number of points Q_(k) used for theB-spline interpolation while still providing good accuracy in thereconstruction. Optionally, user input can be used to define a settingfor the finite difference interval used in the scalar second derivativecomputations. Thus, the finite difference interval size in the scalarsecond derivative calculation can serve as an adjustable tuningparameter which can be controlled to define a quality and efficiency forthe contour reconstruction.

According to a third technique of point reduction for step 604, theshape-salient points are determined as a function of the DeBoor equalenergy theorem. (See DeBoor, C., A Practical Guide to Splines, Springer,New York, 2001, the entire disclosure of which is incorporated herein byreference). With the DeBoor equal energy theorem, the total curvature ofthe entire curve is divided into s equal parts, and the sampled pointsare placed along the curve, at s non-uniform intervals, but in such away as to divide the total curvature into equal parts.

The DeBoor theorem then measures the curvature as the k-th root ofabsolute value of the k-th derivative of the curve,

$\begin{matrix}{{{D^{k}{C(u)}}}^{\frac{1}{k}} = {{{\frac{^{k}}{u^{k}}{C(u)}}}^{\frac{1}{k}}.}} & (13)\end{matrix}$

where D^(k)C(u) denotes the derivative operator. The above-referencedwork by DeBoor proves two instances of a theorem (Theorem II(20),Theorem XII(34)) that optimally places breakpoints (sample points) tointerpolate a curve with minimum error. For a closed curve C(u) oflength L such that 0≦u<L, one can define a set of arc length valuesν_(j)j=1, . . . ,s such that points on C at those values evenly dividethe total curvature. The total curvature K is

$\begin{matrix}{K = {\int_{0}^{L}{{{D^{k}{C(u)}}}^{\frac{1}{k}}{u}}}} & (14)\end{matrix}$

so that dividing it into s equal parts where the energy of any part is1/s of the energy of the curve, or

$\begin{matrix}{{\int_{\upsilon_{j}}^{\upsilon_{j + 1}}{{{D^{k}{C(u)}}}^{\frac{1}{k}}{u}}} = {\frac{1}{s}{\int_{0}^{L}{{{D^{k}{C(u)}}}^{\frac{1}{k}}{{u}.}}}}} & (15)\end{matrix}$

This measure is similar to the ∫(D^(k)C(u))²du “bending energy”curvature measure (see Wahba, G., Spline Models for Observational Data,SIAM (Society for Industrial and Applied Mathematics), Philadelphia,Pa., 1990, the entire disclosure of which is incorporated herein byreference) minimized by spline functions, so it is deemed appropriate tocall the DeBoor technique described herein as an “equal energy” theoremor method.

FIG. 11 shows the plot 1102 of the DeBoor energy versus arc length, anda plot 1104 of the cumulative distribution of that energy for the randomshape curve 400 of FIG. 4. In this plot, the finite differences weretaken across 10% of the curve length.

FIG. 12 shows three DeBoor energy reconstructions of the random shapecurve 400 of FIG. 4 (shown as dotted lines in FIG. 12). The number ofpoints Q_(k) for each reconstruction is set by the chosen value for s.For larger values of s, larger values of Q_(k) will result. Forreconstructed contour 1202, the number of points Q_(k) was 7. Forreconstructed contour 1204, the number of points Q_(k) was 10. Forreconstructed contour 1206, the number of points Q_(k) was 20. As can beseen from FIG. 12, the DeBoor equal energy function producesincreasingly accurate reconstructions of the random figures as thenumber of points s is increased. Unlike the curvature and scalar secondderivative methods described above, the DeBoor equal energy method fitsa given number of points to the curve, instead of the number of pointsbeing set by the detection of underlying peaks in the curve. It shouldbe noted, however, for shapes with a high shape complexity or a largevalue of the integral (14), it may be difficult to anticipate theappropriate number of Q_(k). Optionally, user input can be used todefine a setting for s used in the DeBoor computations. Thus, theinterval size s in the DeBoor equal energy method can serve as anadjustable tuning parameter which can be controlled to define a qualityand efficiency for the contour reconstruction.

Returning to FIG. 6, after a reduced set of points has been generated atstep 604 for each initial set of input points 500, step 606 operates tore-order the points within each reduced set with respect to a frame ofreference and a direction of rotation. Preferably, the frame ofreference is set so that each reduced set's origin (i.e., first) pointis at the cranial-most position on the contour. Also, the direction ofrotation is preferably set to be clockwise, as determined by the methodof Turning Tangents. (See the above-referenced work by DoCarmo). In thismanner, for a given reduced point set, the origin point is set as thecranial-most point, and subsequent points are ordered based on aclockwise rotation starting from the origin point. However, it should beunderstood that other frames of reference and/or directions of rotationcould be used.

Next, step 608 operates to find the points of intersection within eachreduced set of re-ordered points on a desired new plane (e.g., a planethat is non-parallel to the at least one viewing plane for the contoursdefined by the reduced sets of input points). Preferably, step 608operates to find the points within each reduced re-ordered point setthat intersect a desired T plane.

After the points of intersection in the new plane (e.g., a T plane) arefound, step 610 operates to generate a new contour in this new plane byordering the points of intersection and interpolating through the pointsof intersection using B-spline interpolation as described above inSection III.

Thereafter, at step 612, a comparison can be made between the newcontour generated at step 610 and a corresponding patient image in thesame plane. Such a comparison can be made visually by a user. If thegenerated contour is deemed a “match” to the image (i.e., a closecorrespondence between the generated contour and the correspondinganatomy in the displayed image), then the generated contour can bearchived for later use (step 614). If the generated contour is notdeemed a match to the image, then process flow of FIG. 6 can begin anewwith new input points, wherein these new input points perhaps definemore contours than were previously used.

FIGS. 13( a)-(c) depict how a set of B-spline T contours can begenerated from input points corresponding to three input S contours.FIG. 13( a) depicts a GUI 1300 through which a user can define contourswithin various CT image slices of a patient. Different frames of the GUI1300 correspond to different viewing planes into the patient's CT data.Frame 1302 depicts a slice of image data in the T viewing plane. Frame1304 depicts a slice of image data in the C viewing plane, and frame1306 depicts a slice of image data in the S viewing plane. The user cannavigate from slice-to-slice within the GUI 1300 using conventionalsoftware tools. Because the process flow of FIG. 6 allows the user todraw original contours for an anatomical region of interest in anyviewing plane, the user can select the viewing plane(s) in which to drawthe contours based on which viewing plane(s) most clearly depict theanatomical region of interest. In this example, the user has drawn threecontours 1308, 1310, and 1312 in the S viewing plane. The correspondingfootprints for these three S contours are shown in frames 1302 and 1304for the T and C viewing planes respectively. Also, the upper right handframe of GUI 1300 depicts perspective views of these three S contours.The process flow of FIG. 6 can be invoked to generate a plurality of Tcontours 1320 from the sampled points for the three S contours 1308,1310, 1312. FIG. 13( b) depicts several of these generated T contours1320. Each generated T contour 1320 corresponds to a T contour generatedfrom B-spline interpolation as described in connection with FIG. 6 for adifferent T plane.

Furthermore, as can be seen in FIG. 13( b), the generated T contourswrap the three original S contours, and a 3D surface 1330 can berendered from these S and T contours using a series of B-splineinterpolations as described above. FIG. 13( c) depicts a rotated view ofthe contours and surface rendering from FIG. 13( b).

In the embodiment of FIG. 5( b), the input points 500 are representativeof a plurality of any combination of T, S, and/or C contours.Furthermore, as described hereinafter, the software program executed byprocessor 502 in the FIG. 5( b) embodiment is configured to generate a3D surface 506 from the input points 500. From this 3D surface 506,contours of any obliquity, including T contours, can be readilygenerated.

FIG. 14 depicts an exemplary process flow for the FIG. 5( b) embodiment.At step 1402, the software program receives the plurality of inputpoints 500 as described above in connection with step 602 of FIG. 6. Aswith step 602, these input points 500 can be grouped as a plurality ofdifferent initial sets of data points 500, wherein each initial pointset is representative of a different contour. It should be noted thatfor a preferred embodiment of the process flow of FIG. 6, the FIG. 6process flow generates T contours from input contours in the S and/or Cviewing planes. However, with a preferred embodiment of the process flowof FIG. 14, the input contours can be in any viewing plane, including Tcontours.

Next, at step 1404, each initial point set is processed to generate areduced set of input points, as described above in connection with step604 of FIG. 6. By compressing the number of points used to represent thedifferent input contours, the computation of the variational implicitsurface becomes practical. Without such compression of the contourrepresentations, the computational time for generating the variationalimplicit surface from the contour representations requires an undueamount of time on conventional computing resources.

Thereafter, at step 1406, a variational implicit surface is generatedfrom the reduced sets of data points. The variational implicit surfaceis a solution to the scattered data interpolation problem in which thegoal is to determine a smooth function that passes through discrete datapoints. (See Turk and O'Brien, Shape Transformation Using VariationalImplicit Functions, Proceedings of SIGGRAPH 99, Annual ConferenceSeries, pp. 335-342, Los Angeles, Calif., August 1999, the entiredisclosure of which is incorporated herein by reference). For a set ofconstraint points {c₁, . . . ,c_(k)} with a scalar height {h₁, . . .,h_(k)} at each position, one can determine a functionf(x),x=(x,y,z)^(T) that passes through each c_(i) such thatf(c_(i))=h_(i). A variational solution that minimizes the so-called“bending energy” (see the above-referenced work by Turk and O'Brien) isthe sum

$\begin{matrix}{{f(x)} = {{\sum\limits_{j = 1}^{n}{d_{j}{\varphi \left( {x - c_{j}} \right)}}} + {P(x)}}} & (16)\end{matrix}$

over radial basis functions φ_(j) (described below) weighted by scalarcoefficients d_(j), and where c_(j) are the constraint point locationsand P(x) is a degree one polynomial

P(x)=p ₀ +p ₁ x+p ₂ y  (17)

that accounts for constant and linear parts of the function f(x). Theradial basis functions for the 3D constraints appropriate for thisproblem are

φ(x)=|x³|.  (18)

Solving for the constraints h_(i) in terms of the known positions

$\begin{matrix}{h_{i} = {{\sum\limits_{j = 1}^{n}{d_{j}{\varphi \left( {c_{i} - c_{j}} \right)}}} + {P\left( c_{i} \right)}}} & (19)\end{matrix}$

gives a linear system that for 3D constraints c_(i)=(c_(i) ^(x), c_(i)^(y), c_(i) ^(z)) is

$\begin{matrix}{{\begin{bmatrix}\varphi_{11} & \varphi_{12} & \ldots & \varphi_{1k} & 1 & c_{1}^{x} & c_{1}^{y} & c_{1}^{z} \\\varphi_{21} & \varphi_{22} & \ldots & \varphi_{2k} & 1 & c_{2}^{x} & c_{2}^{y} & c_{2}^{z} \\\vdots & \vdots & \; & \vdots & 1 & \vdots & \vdots & \vdots \\\varphi_{k\; 1} & \varphi_{k\; 2} & \ldots & \varphi_{kk} & 1 & c_{k}^{x} & c_{k}^{y} & c_{k}^{z} \\1 & 1 & \ldots & 1 & 0 & 0 & 0 & 0 \\c_{1}^{x} & c_{2}^{x} & \ldots & c_{k}^{x} & 0 & 0 & 0 & 0 \\c_{1}^{y} & c_{2}^{y} & \ldots & c_{k}^{y} & 0 & 0 & 0 & 0 \\c_{1}^{z} & c_{2}^{z} & \ldots & c_{k}^{z} & 0 & 0 & 0 & 0\end{bmatrix}\begin{bmatrix}d_{1} \\d_{2} \\\vdots \\d_{k} \\p_{0} \\p_{1} \\p_{2} \\p_{3}\end{bmatrix}} = \begin{bmatrix}h_{1} \\h_{2} \\\vdots \\h_{k} \\0 \\0 \\0 \\0\end{bmatrix}} & (20)\end{matrix}$

This system is symmetric and positive semi-definite, so there willalways be a unique solution for the d_(j) and the p₁. The solution canbe obtained using LU decomposition. (See the above-referenced works byPress et al. and Golub and Van Loan). In a preferred embodiment, theimplementation of LU decomposition can be the LAPACK implementation thatis known in the art. (See Anderson et al., LAPACK User's Guide, ThirdEdition, SIAM—Society for Industrial and Applied Mathematics,Philadelphia, 1999, the entire disclosure of which is incorporatedherein by reference).

A further feature of the variational implicit surface computation asdescribed in the above-referenced work by Turk and O'Brien is the use ofadditional constraint points, located off the boundary along normalsconnecting with the on-boundary constraints, to more accurately andreliably interpolate the surface through the on-boundary constraints. Ina preferred embodiment, the on-boundary constraints' h_(j) values can beset to 0.0 and the off-boundary values can be set to 1.0. However, asshould be understood, other values can be used in the practice of thisembodiment of the invention.

FIG. 15 demonstrates, for a pair of actual manually drawn contours 1502and 1504, the non-uniform data sampling, the results of the DeBoorenergy analysis using 20 points with both contours, and the constructionof a complete set of constraints for input to the implicit functioncomputation. FIG. 15 depicts the original contour points 1506 for an Scontour 1502 and a C contour 1504. Also depicted in FIG. 15 are theshape-salient points 1508 (shown as the darker points along thecontours) computed from the original points 1506 using theabove-described DeBoor equal energy technique. Furthermore, FIG. 15depicts the normals 1510 (shown as boxes) to the on-contour constraintpoints 1508. The implicit function constraint points would thus includeboth the on-contour shape-salient points 1508 and their correspondingnormals 1510.

Performance of this solution depends partly on the form of the radialbasis function φ(x)one uses, and on the size of the system parameter k(number of all constraint points). The performance of the LU solution ofequation (20) can be done using different choices of Φ. (See Dinh, etal., Reconstructing surfaces by volumetric regularization using radialbasis functions, IEEE Transactions on Pattern Analysis and MachineIntelligence, 24, pp. 1358-1371, 2002, the entire disclosure of which isincorporated herein by reference).

The function |x³| is monotonic increasing, meaning that the matrix in(20) has large off-diagonal values for all constraint point pairs c_(i),c_(j), i≠j. To make the linear system perform more robustly, theabove-referenced work by Dinh describes a modification of the system tomake it more diagonally dominant by adding to the diagonal elements aset of scalar values λ_(i)

$\begin{matrix}{{\begin{bmatrix}{\varphi_{11} + \lambda_{1}} & \varphi_{12} & \ldots & \varphi_{1k} & 1 & c_{1}^{x} & c_{1}^{y} & c_{1}^{z} \\\varphi_{21} & {\varphi_{22} + \lambda_{2}} & \ldots & \varphi_{2k} & 1 & c_{2}^{x} & c_{2}^{y} & c_{2}^{z} \\\vdots & \vdots & \; & \vdots & 1 & \vdots & \vdots & \vdots \\\varphi_{k\; 1} & \varphi_{k\; 2} & \ldots & {\varphi_{kk} + \lambda_{k}} & 1 & c_{k}^{x} & c_{k}^{y} & c_{k}^{z} \\1 & 1 & \ldots & 1 & 0 & 0 & 0 & 0 \\c_{1}^{x} & c_{2}^{x} & \ldots & c_{k}^{x} & 0 & 0 & 0 & 0 \\c_{1}^{y} & c_{2}^{y} & \ldots & c_{k}^{y} & 0 & 0 & 0 & 0 \\c_{1}^{z} & c_{2}^{z} & \ldots & c_{k}^{z} & 0 & 0 & 0 & 0\end{bmatrix}\begin{bmatrix}d_{1} \\d_{2} \\\vdots \\d_{k} \\p_{0} \\p_{1} \\p_{2} \\p_{3}\end{bmatrix}} = {\begin{bmatrix}h_{1} \\h_{2} \\\vdots \\h_{k} \\0 \\0 \\0 \\0\end{bmatrix}.}} & (20)\end{matrix}$

A preferred embodiment uses values λ_(Boundary)=0.001 andλ_(OffBoundary)=1.0. However, it should be understood that other valuescould be used.

After solving for the d_(j) and the p_(j) in Equation (20), the implicitfunction in (16) can be evaluated to determine that set of points{x_(i)} for which f(x_(i))=0. (The zero-th level of f(x) is that onwhich the boundary points lie). The method of Bloomenthal (seeBloomenthal, J., An Implicit Surface Polygonizer, Graphics Gems IV, P.Heckbert, Ed., Academic Press, New York, 1994, the entire disclosure ofwhich is incorporated herein by reference) can be used to track aroundthe function and determine the locations of mesh nodes from which a 3Dsurface may be constructed. A closed surface constructed in this way canbe termed a variational implicit surface. (See the above-referenced workby Turk and O'Brien).

At step 1408, that mesh can then be clipped by planes parallel to thexz-plane at the appropriate y-value(s) to produce the desired Tcontour(s) for display to the user. The mesh representation and clippingfunctionality can be performed using the VTK software system availablefrom Kitware, Inc. of Clifton Park, N.Y. (See Schroeder et al., TheVisualization Toolkit, 4^(th) Ed., Kitware, 2006, the entire disclosureof which is incorporated herein by reference).

Thereafter, as with steps 612 and 614 of FIG. 6, the generated contourcan be compared to the image (step 1410) and archived if it sufficientlymatches the anatomy of interest shown in the image (step 1412). If not asufficient match, the process of FIG. 14 can begin anew.

As indicated, the main performance limitation for computing avariational implicit surface is the total number of constraints, and fork greater than a few thousand, the variational implicit surfacecomputation takes too much time to be useful for real-time applications.However, the inventor herein believes that by reducing the number ofconstraint points used for computing the variational implicit surfacevia any of the compression operations described in connection with steps1404 and 604 for contour representations, the computation of variationalimplicit surfaces will become practical for 3D medical contouring.Furthermore, a fortunate property of the variational implicit surface isits ability to forgive small mismatches in orthogonal contours that arerequired to intersect (because they are curves on the same surface) butdo not because the user was unable to draw them carefully enough. Forexample, when the sampling interval in the Bloomenthal algorithm is setto the inter-T plane distance, the resulting surfaces are sampled at toocoarse a level to reveal the small wrinkles in the actual surface, andthe resulting T contours are not affected by the missed T/S/Cintersections.

FIG. 16 illustrates the variational implicit surface 1602 produced fromthe contours 1502 and 1504 in FIG. 15. FIG. 16 also depicts cutaways1604 and 1606 that show the profiles of the contours 1502 and 1504,respectively, that were used to compute surface 1602.

FIGS. 17-19( b) demonstrate the application of the variational implicitsurface method in the contouring of three organs shown in FIG. 1—theprostate 110, bladder 112, and rectum 114. In FIG. 17, variationalimplicit contours of the prostate 110 are shown in wireframe along witha conventional T-only contoured rendering. The top part of FIG. 17depicts a left sagittal view of the variational implicit contours. Thebottom part of FIG. 17 depicts a frontal (anterior) coronal view of thevariational implicit contours. The wireframe for FIG. 17 was createdusing a single C contour and three S contours as inputs.

In FIG. 18, variational implicit contours of the bladder 112 are shownin wireframe along with a conventional T-only contoured rendering. Thetop portion of FIG. 18 depicts a right sagittal view of the variationalimplicit contours, and the bottom portion of FIG. 18 depicts a frontal(anterior) coronal view of the variational implicit contours. Thewireframe for FIG. 18 was generated using a single C contour and three Scontours as inputs.

In FIGS. 19( a) and (b), variational implicit contours of the rectum 114are shown in wireframe along with a conventional T-only contouredrendering. FIG. 19( a) depicts a left sagittal view of the variationalimplicit contours, and FIG. 19( b) depicts a frontal (anterior) coronalview of the variational implicit contours. The wireframe for FIGS. 19(a) and (b) was generated using a three C contours and five S contours asinputs.

As shown by FIGS. 17-19( b), the wireframes show good agreement with theconventionally drawn structures depicted at the right in these figures.

It should also be noted that it may sometimes be the case wherein aninitial set of input points corresponding to a contour contains only asmall number of points, long gaps in the sequence of input points,and/or two or more points having the same coordinates (x,y). Forexample, both a small number of points and long gaps between pointswould likely result when a user defines a contour by only picking pointsat the vertices of a polygon that approximates the contour. Duplicatepoints can result when the user picks points along a contour because agraphics subsystem will sometimes interpret a single mouse button pushas multiple events. In such instances, the process flow of FIG. 20 canbe employed. At step 2002, an input point set for a contour is received.At step 2004, a check is made as to whether there are points withduplicate (x,y) coordinates. If yes at step 2004, the program proceedsto step 2006, where all duplicate points are removed from the contour'spoint set. If not, the input point set is retained and the processproceeds to step 2008. At step 2008, the program determines whetherthere are gaps in the input point set greater than a threshold fractionof the total contour length. If yes at step 2008, the program proceedsto step 2010 which augments the input point set by filling the gaps withpoints by linear interpolation between the input points at the endpoints of the gap. Optionally, the identities of the original pointsfrom the input point set and the identifies of the points added at step2010 can be preserved so that the process of finding the shape-salientpoints within the point set will only allow the original points to beclassified as shape-salient. If no at step 2008, then the programproceeds to step 2012. At step 2012, a contour is interpolated from thepoints in the input point set using B-spline interpolation as previouslydescribed in Section III. Thereafter, at step 2014, a plurality ofpoints can be sampled from the interpolated contour, and these sampledpoints can be used to replace and/or augment the points in the inputpoint set used to represent that contour. Step 2016 then operates tocheck whether the user has defined any additional contours. If not, theprocess flow can proceed to further complete the contouring operations(such as by proceeding to step 604 of FIG. 6 or step 1404 of FIG. 14).If the user has provided further input, the process can return to step2002.

It should also be noted that the B-spline interpolation and variationalimplicit surface generation can be combined in a single process flow asdifferent modes of operation, as shown in FIG. 21. At step 2102, variousinput point sets for different contours are received. Then, at step2104, the process flow decides which reconstruction mode should beused—e.g., a B-spline interpolation mode or a variational implicitsurface mode. The decision at step 2104 can be made in any of a numberof ways. For example, the B-spline interpolation mode can be used whereonly a small number (e.g., less than or equal to three) of S and/or Cinput contours have been defined, and the variational implicit surfacemode can be used for other cases. Further still, a user can select whichmode to be used via some form of user mode input.

If the B-spline interpolation mode is used, then steps 2106 and 2108 canbe performed, wherein these steps correspond to steps 606 and 608 fromFIG. 6, albeit without the preceding point reduction operation of step604. However, it should also be noted that steps 2106 and 2108 could bereplaced by steps 604, 606, and 608 if desired by a practitioner of thisembodiment of the invention. Step 610 preferably operates as describedabove in connection with FIG. 6.

If the variational implicit surface mode is used, then steps 1404, 1406,and 1408 can be followed as described in connection with FIG. 14.

While the present invention has been described above in relation to itspreferred embodiments, various modifications may be made thereto thatstill fall within the invention's scope. Such modifications to theinvention will be recognizable upon review of the teachings herein.Accordingly, the full scope of the present invention is to be definedsolely by the appended claims and their legal equivalents.

What is claimed is:
 1. A computer-implemented method for generating athree-dimensional (3D) surface corresponding to a region of interestwithin an image, the method comprising: generating a second plurality ofdata points from a first plurality of data points, the second pluralitybeing less than the first plurality, the first plurality of data pointsbeing representative of a plurality of contours corresponding to theregion of interest; and generating a 3D surface representative of theregion of interest based on the second plurality of data points; andwherein the method steps are performed by a processor.
 2. The method ofclaim 1 wherein the 3D surface generating step comprises the processorgenerating a variational implicit surface representative of the regionof interest based on the second plurality of data points.
 3. The methodof claim 2 wherein the step of generating the second plurality of pointscomprises: the processor reducing the first plurality of points to aplurality of shape-salient points.
 4. The method of claim 2 wherein thestep of generating the second plurality of data points comprises: theprocessor computing a plurality of curvature values for the firstplurality of data points; and the processor defining the secondplurality of data points as a function of the computed curvature values.5. The method of claim 4 wherein the defining step comprises theprocessor choosing as the second plurality of data points those pointswithin first plurality of data points for which the computed curvaturevalue satisfies a predetermined condition.
 6. The method of claim 2wherein the step of generating the second plurality of data pointscomprises: the processor computing a plurality of scalar secondderivative values for the first plurality of data points; and theprocessor defining the second plurality of data points as a function ofthe computed scalar second derivative values.
 7. The method of claim 6wherein the defining step comprises the processor choosing as the secondplurality of data points those points within first plurality of datapoints for which the computed second scalar derivative value satisfies apredetermined condition.
 8. The method of claim 2 wherein the step ofgenerating the second plurality of data points comprises: the processordefining the second plurality of data points based on a DeBoor equalenergy theorem function.
 9. The method of claim 8 wherein the definingstep comprises: the processor computing a total curvature for eachcontour defined by the first plurality of data points; the processordividing the total curvature for each contour defined by the firstplurality of data points into a plurality s of equal energy parts; andthe processor selecting as the second plurality of data points the endpoints of each of the s arc lengths that result from the dividing step.10. The method of claim 1 further comprising: the processor controllinghow many of the second plurality of data points are generated via auser-adjustable tuning parameter.
 11. An apparatus for generating athree-dimensional (3D) surface corresponding to a region of interestwithin an image, the apparatus comprising: a processor configured to (1)generate a second plurality of data points from a first plurality ofdata points, the second plurality being less than the first plurality,the first plurality of data points being representative of a pluralityof contours corresponding to the region of interest, and (2) generate a3D surface representative of the region of interest based on the secondplurality of data points.
 12. The apparatus of claim 11 wherein theprocessor is further configured to generate the 3D surface as avariational implicit surface.
 13. The apparatus of claim 12 wherein theprocessor is further configured to generate the second plurality ofpoints by reducing the first plurality of points to a plurality ofshape-salient points.
 14. The apparatus of claim 12 wherein theprocessor is further configured to (1) compute a plurality of curvaturevalues for the first plurality of input points, and (2) define thesecond plurality of data points as a function of the computed curvaturevalues.
 15. The apparatus of claim 12 wherein the processor is furtherconfigured to (1) compute a plurality of scalar second derivative valuesfor the first plurality of input points, and (2) define the secondplurality of data points as a function of the computed scalar secondderivative values.
 16. The apparatus of claim 12 wherein the processoris further configured to define the second plurality of data pointsbased on a DeBoor equal energy theorem function.
 17. A computer-readablestorage medium for generating a three-dimensional (3D) surfacecorresponding to a region of interest within an image, thecomputer-readable storage medium comprising: a plurality ofcomputer-executable instructions for (1) generating a second pluralityof data points from a first plurality of data points, the secondplurality being less than the first plurality, the first plurality ofdata points being representative of a plurality of contourscorresponding to the region of interest, and (2) generating a 3D surfacerepresentative of the region of interest based on the second pluralityof data points, and wherein the instructions are resident on thecomputer-readable storage medium.
 18. The computer-readable storagemedium of claim 17 wherein the instructions further comprisecomputer-executable instructions for generating the 3D surface as avariational implicit surface.
 19. The computer-readable storage mediumof claim 18 wherein the instructions further comprisecomputer-executable instructions for generating the second plurality ofpoints by reducing the first plurality of points to a plurality ofshape-salient points.
 20. The computer-readable storage medium of claim18 wherein the instructions further comprise computer-executableinstructions for (1) computing a plurality of curvature values for thefirst plurality of input points, and (2) defining the second pluralityof data points as a function of the computed curvature values.
 21. Thecomputer-readable storage medium of claim 18 wherein the instructionsfurther comprise computer-executable instructions for (1) computing aplurality of scalar second derivative values for the first plurality ofinput points, and (2) defining the second plurality of data points as afunction of the computed scalar second derivative values.
 22. Thecomputer-readable storage medium of claim 18 wherein the instructionsfurther comprise computer-executable instructions for defining thesecond plurality of data points based on a DeBoor equal energy theoremfunction.
 23. A computer-implemented method for processing a pluralityof data point sets, each data point set comprising a plurality of datapoints that are representative of a contour corresponding to a region ofinterest within an image volume, each contour having a correspondingplane, the method comprising: for at least one of the data points sets,interpolating a contour from the points within the at least one datapoint set; determining a plurality of shape-salient points from theinterpolated contour; and processing a portion of the data points withinthe data point sets to thereby generate a new contour in a new plane,the new contour corresponding to the region of interest, and wherein theportion includes the determined shape-salient points; and wherein themethod steps are performed by a processor.
 24. A computer-implementedmethod for processing contour data, the method comprising: performing acompression operation on a plurality of contour representations tothereby generate a plurality of compressed contour representations, eachcontour representation comprising data corresponding to a contour of aregion of interest within an image, each contour residing in a plane;and processing the compressed contour representations to therebygenerate a new contour in a new plane, the new contour corresponding tothe region of interest; and wherein the method steps are performed by aprocessor.
 25. A computer-implemented method for processing a pluralityof data point sets, each data point set comprising a plurality of datapoints that are representative of a contour corresponding to a region ofinterest within an image volume, each contour having a correspondingplane, the method comprising: for at least one of the data points sets,interpolating a contour from the points within the at least one datapoint set; determining a plurality of shape-salient points from theinterpolated contour; and processing a portion of the data points withinthe data point sets to thereby generate a variational implicit surfacethat corresponds to the region of interest, and wherein the portionincludes the determined shape-salient points; and wherein the methodsteps are performed by a processor.
 26. The method of claim 25 whereinthe interpolating step comprises the processor interpolating the contourfrom the points within the at least one data point set using B-splineinterpolation.
 27. A computer-implemented method for processing contourdata, the method comprising: performing a compression operation on aplurality of contour representations to thereby generate a plurality ofcompressed contour representations, each contour representationcomprising data corresponding to a contour of a region of interestwithin an image, each contour residing in a plane; and processing thecompressed contour representations to thereby generate a variationalimplicit surface corresponding to the region of interest; and whereinthe method steps are performed by a processor.